VIB/IR analysis workflow using MatSciToolkit

VIB/IR analysis workflow using MatSciToolkit#

This tutorial will guide you through the use of the MatSciToolkit package to calculate the vibrational frequencies of a molecule or surface. The MatSciToolkit package is a python package the provides a set of tools to automate the workflow for the vibrational and IR analysis. The package is built on top of the Atomic Simulation Environment (ASE) and simplifies the workflow while using the Quantum Espresso code as the calculator.

A quick example of the data we get is the following

Hide code cell content
from ase.build import molecule
from ase.calculators.emt import EMT
from ase.optimize import BFGS
from ase.vibrations import Vibrations
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
import matplotlib.animation as animation
h2o = molecule('H2O', cell=(10,10,10))
h2o.center()
h2o.calc = EMT()
BFGS(h2o).run(fmax=0.01)
Hide code cell output
      Step     Time          Energy         fmax
BFGS:    0 16:24:28        2.619811        7.7384
BFGS:    1 16:24:28        1.912906        1.3454
BFGS:    2 16:24:28        1.882033        0.4035
BFGS:    3 16:24:28        1.879275        0.0317
BFGS:    4 16:24:28        1.879253        0.0096
True
Hide code cell source
import plotly.graph_objects as go
from ase.data.colors import jmol_colors
from ase.data import atomic_numbers, covalent_radii
import numpy as np

def ms(x, y, z, radius, resolution=20):
    """Return the coordinates for plotting a sphere centered at (x,y,z)"""
    u, v = np.mgrid[0:2*np.pi:resolution*2j, 0:np.pi:resolution*1j]
    X = radius * np.cos(u)*np.sin(v) + x
    Y = radius * np.sin(u)*np.sin(v) + y
    Z = radius * np.cos(v) + z
    return (X, Y, Z)


pos = h2o.get_positions()
numbers = h2o.get_atomic_numbers()
colors = [jmol_colors[number] for number in numbers]
cs = np.array(colors) * 255
cs = [f"rgb({int(color[0])}, {int(color[1])}, {int(color[2])})" for color in cs]
scale = 1
radii = covalent_radii[numbers]
radii = [(scale * r)**1 for r in radii]

data = []
for i in range(len(pos)):
    (x_pns_surface, y_pns_surface, z_pns_surface) = ms(pos[i,0], pos[i,1], pos[i,2], radii[i])
    data.append(go.Surface(x=x_pns_surface, y=y_pns_surface, z=z_pns_surface, opacity=1, showscale=False, colorscale=[cs[i], cs[i]]))

fig = go.Figure(data=data)

fig.update_layout(
    title='Water molecule structure',
    scene=dict(
        xaxis=dict(title='X', range=[0, 10]),
        yaxis=dict(title='Y', range=[0, 10]),
        zaxis=dict(title='Z', range=[0, 10]),
        aspectmode='cube'
    ),
    width=1000,  # Set the width of the figure
    height=500,  # Set the height of the figure
    margin=dict(l=0, r=0, t=30, b=0),  # Set margin to 0 for a tight layout
    template='plotly_dark',
    
)

fig.show()

We’ll get the vibrations or infrared data#

vib = Vibrations(h2o)
vib.run()
vib.summary()
---------------------
  #    meV     cm^-1
---------------------
  0    6.3i     51.0i
  1    5.9i     47.6i
  2    0.0i      0.3i
  3    0.0i      0.1i
  4    0.1       0.4
  5    5.4      43.8
  6   32.1     258.9
  7  296.7    2392.7
  8  387.4    3124.9
---------------------
Zero-point energy: 0.361 eV
Hide code cell source
vibdata = vib.get_vibrations()

To visualize, the mode 6, 7, and 8 is:#

Hide code cell source
%matplotlib notebook
from IPython.display import HTML
plt.ioff()
fig, (ax1, ax2, ax3) = plt.subplots(3, 2, figsize=(8, 8), dpi=150)
axs = [ax1, ax2, ax3]

itervib6 = list(vibdata.iter_animated_mode(6))
itervib7 = list(vibdata.iter_animated_mode(7))
itervib8 = list(vibdata.iter_animated_mode(8))
itervibs = [itervib6, itervib7, itervib8]

fig.tight_layout()
def animate(i):
    for ax in axs:
        ax[0].cla()
        ax[1].cla()
        
    for k, (ax, itervib) in enumerate(zip(axs, itervibs)):
        plot_atoms(itervib[i], ax=ax[0], radii=0.9, rotation=(''))
        plot_atoms(itervib[i], ax=ax[1], radii=0.9, rotation=('90y'))
        
        mode = {0: "Mode 6", 1: "Mode 7", 2: "Mode 8"}
        ax[0].set_title(mode[k])
        ax[1].set_title(mode[k])
    
        # if k == 2: 
        ax[0].set_xlabel("X-axis")
        ax[1].set_xlabel("Z-axis")
    
        ax[0].set_ylabel("Y-axis")
        ax[1].set_ylabel("Y-axis")

    if i == 0:
        fig.tight_layout()
        
ani = animation.FuncAnimation(fig, animate, frames=len(itervib6), interval=25)


# from IPython.display import HTML, Javascript
HTML(ani.to_jshtml())
# HTML(ani.to_html5_video())